% Table 4 and Figure 1 in "Commodity Prices, Convenience Yields, and
% Inflation," March 2011
% Note: The estimation of IMA(1,1) uses the GARCH toolbox.

clear all;

load cyj.dat;                   % convenience yields on 23 commodities
load qj_onesided.dat;           % 1-sided detrended real commodity prices
load us_cpi_monthly.dat;        % US CPI indeces (all items, lfe etc.)

qj=qj_onesided;
P=us_cpi_monthly;
nn=rows(P);

rhat=2;         % number of principal components
nwl=4;          % number of lags in Newey-West HAC estimator
maxp=4;         % maximum lag for AIC
fobs=24;        % # of observations for one-sided MA
hh=[1 3 6 12];  % forecast horizon

cyj=standard(cyj);
cyj=cyj(fobs:nn,:);
P=P(fobs:nn,:);
nn=rows(cyj);

for inflmeas=1:2;     % inflation measure: (1) all items, (2) lfe

for hi=1:4;
h=hh(hi);

Spec = garchset('R', 0, 'M', 1, 'Display', 'off');

TT=154;
act=[]; forc1=zeros(nn-TT-h,maxp); forc2=zeros(nn-TT-h,maxp); 
forc3=zeros(nn-TT-h,maxp); forc4=[]; forc5=[]; 
mse1=[]; mse2=[]; mse3=[]; mse4=[]; mse5=[]; 
mae1=[]; mae2=[]; mae3=[]; mae4=[]; mae5=[]; 
for j=1:nn-TT-h;
    infl_h=100*(log(P(1+maxp+h:TT+j,inflmeas))-log(P(1+maxp:TT+j-h,inflmeas)));
    dinfl = 100*(log(P(1+maxp:TT+j,inflmeas))-log(P(maxp:TT+j-1,inflmeas)))-...
            100*(log(P(maxp:TT+j-1,inflmeas))-log(P(maxp-1:TT+j-2,inflmeas)));
    [ehat,Fhat,lamhat,ve]=pc(standard(cyj(1:TT+j,:)),rhat);
    [ehat2,Fhat2,lamhat2,ve2]=pc(standard(qj(1:TT+j,:)),rhat);

    const=ones(rows(infl_h),1);
    nr=rows(const);
    aic1=zeros(maxp,1); aic2=zeros(maxp,1); aic3=zeros(maxp,1);
    infl_le=zeros(nr,maxp); p1_cy=zeros(nr,maxp); p2_cy=zeros(nr,maxp);
    p1_rcpd=zeros(nr,maxp); p2_rcpd=zeros(nr,maxp); roil=zeros(nr,maxp);
    infl_lf=zeros(1,maxp); p1_cyf=zeros(1,maxp); p2_cyf=zeros(1,maxp);
    p1_rcpdf=zeros(1,maxp); p2_rcpdf=zeros(1,maxp); roilf=zeros(1,maxp);
    for ip=1:maxp;
        infl_le(:,ip)=100*(log(P(2+maxp-ip:TT+j-h+1-ip,inflmeas))-log(P(1+maxp-ip:TT+j-h-ip,inflmeas)));
        p1_cy(:,ip)=Fhat(2+maxp-ip:TT+j-h+1-ip,1);
        p2_cy(:,ip)=Fhat(2+maxp-ip:TT+j-h+1-ip,2);
        p1_rcpd(:,ip)=Fhat2(2+maxp-ip:TT+j-h+1-ip,1);
        p2_rcpd(:,ip)=Fhat2(2+maxp-ip:TT+j-h+1-ip,2);
        roil(:,ip)=qj(2+maxp-ip:TT+j-h+1-ip,22);
        infl_lf(:,ip)=100*(log(P(TT+j+1-ip,inflmeas))-log(P(TT+j-ip,inflmeas)));
        p1_cyf(:,ip)=Fhat(TT+j+1-ip,1);
        p2_cyf(:,ip)=Fhat(TT+j+1-ip,2);
        p1_rcpdf(:,ip)=Fhat2(TT+j+1-ip,1);
        p2_rcpdf(:,ip)=Fhat2(TT+j+1-ip,2);
        roilf(:,ip)=qj(TT+j+1-ip,22);
        res=nwest(infl_h,[const p1_cy(:,1:ip) p2_cy(:,1:ip) infl_le(:,1:ip)],nwl);
        res2=nwest(infl_h,[const p1_rcpd(:,1:ip) p2_rcpd(:,1:ip) roil(:,1:ip) infl_le(:,1:ip)],nwl);
        res3=nwest(infl_h,[const infl_le(:,1:ip)],nwl);
        aic1(ip)=log(res.sige)+2*(rows(res.beta)-1)./rows(const);
        aic2(ip)=log(res2.sige)+2*(rows(res2.beta)-1)./rows(const);
        aic3(ip)=log(res3.sige)+2*rows((res3.beta)-1)./rows(const);
        indep1=[1 p1_cyf(:,1:ip) p2_cyf(:,1:ip) infl_lf(:,1:ip)];
        indep2=[1 p1_rcpdf(:,1:ip) p2_rcpdf(:,1:ip) roilf(:,1:ip) infl_lf(:,1:ip)];
        indep3=[1 infl_lf(:,1:ip)];
        forc1(j,ip)=indep1*res.beta;
        forc2(j,ip)=indep2*res2.beta;
        forc3(j,ip)=indep3*res3.beta;
    end;
    [aicmin1,minind1] = min(aic1);
    [aicmin2,minind2] = min(aic2);
    [aicmin3,minind3] = min(aic3);
    res4=nwest(infl_h,const,nwl);
    indep4=[1];
    forc4(j,:)=indep4*res4.beta;
    [Coeff, Errors, LLF, Innovations, Sigmas, Summary]  = garchfit(Spec, dinfl);
    C = garchget(Coeff, 'C'); 
    MA = garchget(Coeff, 'MA');
    Spec = garchset('R', 0, 'M', 1, 'C', C, 'MA', MA, 'Display', 'off');
    forc5(j,:)=C+Innovations(rows(Innovations))*MA + 100*(log(P(TT+j,inflmeas))-log(P(TT+j-1,inflmeas)));
    act(j,:)=100*(log(P(TT+j+h,inflmeas))-log(P(TT+j,inflmeas)));
    mse1(j,:)=(forc1(j,minind1)-act(j))^2;
    mse2(j,:)=(forc2(j,minind2)-act(j))^2;
    mse3(j,:)=(forc3(j,minind3)-act(j))^2;
    mse4(j,:)=(forc4(j)-act(j))^2;
    mse5(j,:)=(h*forc5(j)-act(j))^2;    
    mae1(j,:)=abs(forc1(j,minind1)-act(j));
    mae2(j,:)=abs(forc2(j,minind2)-act(j));
    mae3(j,:)=abs(forc3(j,minind3)-act(j));
    mae4(j,:)=abs(forc4(j)-act(j));
    mae5(j,:)=abs(h*forc5(j)-act(j));
end;

bench1=sqrt(mean(mse3));
bench2=mean(mae3);
if inflmeas==1;
 fprintf('    All Items\n');
else
 fprintf('    Less Food and Energy\n');
end;
fprintf('    h = %6.4f\n',h);
fprintf('    RMSE\n')
disp([sqrt(mean(mse1))./bench1 sqrt(mean(mse2))./bench1 sqrt(mean(mse5))./bench1])
fprintf('    MAE\n')
disp([mean(mae1)./bench2 mean(mae2)./bench2 mean(mae5)./bench2])

if hi==4;
    plot(act,'-b','LineWidth',1)
    hold on
    plot(forc1(:,minind1),'-rs','MarkerSize',4);
    hold off
    set(gca,'XTick',[2 14 26 38 50 62 74 86 98 110])
    set(gca,'XTickLabel',{'1999' '2000' '2001' '2002' '2003' ...
     '2004' '2005' '2006' '2007' '2008'})
    if inflmeas==1;
        h = legend('actual values','forecasts: PCs',2);
    else
        h = legend('actual values','forecasts: PCs',3);
    end
 set(h,'Interpreter','none')
end

end

end